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Abstract. The gravitational wave detectors currently in operation perform 
the analysis of their scientific data jointly. Concerning the search for bursting 
sources, coherent data analysis methods have been shown to be more efficient. 
In the coherent approach, the data collected by the detectors are time-shifted 
and linearly combined so that the signatures received by each detector add up 
constructively (thus improving the resulting signal-to- noise ratio). This operation 
has to be performed over a sky grid (which determines the sky locations to be 
searched). A limitation of those pipelines is their large computing cost. One of the 
available degrees of freedom to reduce the cost is the choice of the sky grid. Ideally, 
the sky sampling scheme should adapt the angular resolution associated with the 
considered gravitational wave detector network. As the geometry of detector 
network is not regular (the detectors are not equally spaced and oriented), the 
angular resolution varies largely depending on the sky location. We propose here 
a procedure which designs sky grids that permit a complete sky coverage with a 
minimum number of vertices and thus adapt the local resolution. 
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1. Introduction 

According to Einstein's theory of General Relativity, the gravitational interaction 
manifests the geometry of space-time curved by matter. The theory also predicts 
the existence of radiative solutions to the space-time dynamics which are called 
gravitational waves. Several detectors (including LIGO in the US, Virgo and GEO 
in Europe) [1] have been built in the last decade to provide the direct observation 
of the gravitational waves. They are mainly long-baseline Michelson-type laser 
interferometers able to sense the weak strain in the arms of the instrument caused 
by the gravitational waves radiated by astrophysical sources. 

While no gravitational wave has been observed so far by these detectors, a first 
detection is plausible with the sensitivities that will be achieved for the up-coming 
data takings (labeled S6/VSR2, to start in 2009). Combining gravitational wave data 
with other types of observations may provide crucial help [2] as they allow better 
background rejection. For instance, a follow-up search for optical counterparts has 
been suggested in [3]. The search for coincident high-energy neutrinos has also been 
proposed in [4]. 

As these combined searches check the occurence of an (optical or neutrino) event 
in time/spatial coincidence with the gravitational wave triggers, they rely on an 
estimate for the location of the gravitational wave source. 

Several data analysis pipelines have been proposed to search for bursts of 
gravitational waves (typically from star collapses or binary mergers), our primary 
focus here. We will be interested in the ones O [6l [71 [8j which perform the coherent 
processing of the data streams collected by each gravitational wave detectors, as they 
have shown to be particularly efficient. 

In the coherent approach, the data collected by the detectors are time-shifted 
and linearly combined so that the signatures received by each detector add up 
constructively if they originate from a given sky location. The transient signal 
is then searched for in the combined data stream. This operation is repeated by 
scanning many sky locations. The result is a likelihood sky map which quantifies the 
"likelihood" of an actual gravitational wave burst source in a particular sky location. 
The computation of likelihood sky maps requires the discretization of the sky sphere. 
A uniform sky sampling (see e.g., [9]) is usually used in practice. 

A fundamental limitation of the coherent pipelines is their large computational 
cost. One of the available degrees of freedom to reduce the cost is the choice of the 
sky grid. On one hand, the computing cost scales linearly with the number of grid 
points or vertices (the same calculation is repeated for each bin of the sky grid with 
different data and coherent mixing coefficients). On the other hand, the grid should 
not be too coarse not to miss any signals. In this paper, we address this trade-off and 
propose a procedure which designs a sky grid that permits a complete sky coverage with 
a minimum number of vertices. 

This procedure determines the location of the grid vertices on the basis of an 
estimate of the local angular resolution of the gravitational wave detector network. It 
thus takes into account both the specific geometry of the detector network and the 
characteristics of the individual detector antenna patterns. Calibration uncertainties 
and other timing errors may be also folded in. 

We give here a brief outline of the method and connect the various parts to the 
corresponding section of the paper. In section[21 we describe the context of this study. 
We describe the characteristics of the signature left by a passing gravitational wave. 
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We spend special attention on the "travel time" (difference in the time of arrival of 
gravitational waves at each detector with respect to a reference one) and its property 
because it is an important ingredient for the sky sampling problem. 

The direction of a source can be equivalently described in terms of the usual 
spherical coordinates or using the set of travel times. In section [3l we produce a 
first sky grid by sampling regularly in the coordinate system associated to the travel 
times. We propose a robust mapping between the discretized travel times and spherical 
angles. 

The sky grid produced in the first step is "over-sampled" as it is based on timing 
information only. In a second step presented in section |4l we use the information of 
the position and orientation of all the detectors and extract the smallest sub-set of 
vertices that ensures complete sky coverage for a given loss in signal-to-noise ratio. 
We build this sub-set by casting the grid size minimization as a set covering problem 
that can be efhciently solved using a greedy procedure based on a dual Lagrangian 
relaxation heuristic. As an illustration, we produce these sky-grids for an idealized 
detector network and for the LIGO- Virgo detector network. Section [5] concludes this 
paper. 

2. Framework 

Let us consider a network of D gravitational wave detectors. The response of detector 
i to an incoming gravitational wave emitted from the direction ((/>, 9) can be written 
as [S]: 

s,{t) = ^[F,icj,,9yhit-nici),e))], t>t, (1) 

where Fi is the complex-valued antenna pattern (x* denotes the complex conjugate 
of x), and h(t) = h+(t) + ihx{t) is the complex-valued gravitational wave signal 
combining the two polarizations + and x. ti is the time of arrival of gravitational 
wave at detector i. In eq. ([1]), time is defined with respect to a reference which we 
arbitrarily set to be the time at detector 1. With this convention, Ti{(j), 9) = {ti — ti) 
denotes the time taken by the wave to travel from detector 1 to detector i. For this 
reason, we refer to Ti{(f>,9) as travel time. It is proportional to the projected distance 
between detectors 1 and i onto the projection direction of the wave, viz. 

r,(0,0) = i(r, -ri)^w(0,0), (2) 
c 

where r^ is the coordinate vector of detector i, w(0, 9) is the unit wave vector and c 
is the speed of light in free space. Here, bold symbols are (column) vectors and x-^ 
denotes the transpose of x. 

2.1. Link between spherical coordinates and travel times 

In this section, we describe the relationship between the spherical coordinates (0, 9) 
and the set of travel times {Ti((/), 0)}i=i,...,D- This has already been examined in [10] 
and this section summarizes some of the results presented in this article. 

We choose to work in the reference frame presented in Fig. [1] where the expressions 
linking spherical coordinates and travel times are simple. Let Ti = du/c be the time 
of flight between detectors Di and Di (here du — ||ri — ri|| is the distance between 
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Figure 1. Reference frame — We choose here to work in a reference frame 
(described for example in llOp in which expressions hnking spherical coordinates 
and travel times are simple. In this graph, Di corresponds to the position of 
detector i. The origin of the network is chosen to be detector Di which also 
define our time reference. The z axis is defined as the line joining Di and D2. 
The X axis is defined to be the intersection between the plane orthogonal to the 
z axis and the plane containing Di, D2 and D3. Finally the y axis is chosen so 
that the x, y and z axes form a right-handed coordinate system. D\S represents 
the direction of the gravitational wave source. 



Di and Di). With these notations, the travel time can be expressed with respect to 
the spherical coordinates (0, 0) in the reference frame by: 



where a2i is the angle between the lines (-D1-D2) and {DiDt) and ijji is the angle 
between the x axis and the line (DiDixy) (with Di^y being the projection of Di onto 
the xy plane). 

The spherical coordinates in this frame can be obtained by the reciprocal 
expressions: 



with i > 3. These equations display well-known features of triangulation: the zenith 
angle 9 can be determined only by two detectors, however the azimuth angle cj) 
requires at least four detectors to be completely determined (including the sign). Thus 
alternatively, the set of travel times may be seen as "coordinates" equivalent to the 
spherical angles. 

As we will describe later in Sec. [Sj the coordinate system generated by the 
collection of travel times is more convenient to produce sky grids adapted to the 
considered detector network. However, this coordinate system is not "standard" . 
This is why we examine its geometry into some more details now. The 2-sphere is 
completely characterized by the two spherical angles. For a number D of detectors. 



Ti (cos 9 cos a2i + sin a2i cos (0 — ?/'«) sin 9) , 



(3) 




(4) 



(5) 
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the number of travel times is D — 1. We need to identify the corresponding 2-sphere 
in the resulting D — 1- dimensional coordinate space. 

From eq. ([3]), we see that the travel time takes values in [—Ti,Ti]. Therefore, 
the admissible values for the coordinates r = [t2, ■ ■ ■ ,tb]'^ are located inside a 
D — 1-rectangular cuboid. But, for a detector network with D > 2, the travel time 
coordinates do not span this cuboid entirely. For instance, it is impossible to have 
T2 = T2 and T3 — T3 unless Di, D2 and D3 lie on a straight line. From the constraints 
imposed by the physics, i.e. the propagation of the wave and the geometry of the 
network, it is possible to determine the admissible surface spanned by the travel time 
coordinates. 

The equation of this two dimensional admissible surface can be obtained from 
Eq. ([3]) [To]. It can be written in the form 

t'^AdT <Bd D = 3 (6) 

= Bd D = A. (7) 

The expressions of the (symmetric) matrix A^i and Bd were obtained in [T(T for _D = 3 
and D = 4. We give them below for completeness. 

For D — the admissible surface is an ellipse with 



A, = 
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cos a23 
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the admissible surface is an ellipsoid. If A4 
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sin q;24 cos Tp4 cot 0:23 ) , 
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T4 sin a24 



cos V'. 



4, 



T3 sin 0:23 
while B4 = T2 siv? a24 sin^ ip^^. 

With the projection of the ellipsoid onto the (r^, Tj) plane, we retrieve the ellipse 
for the corresponding Di, Di and Dj three-detector network. The two points from 
the two shelves of the ellipsoid which get projected in a single one are associated to 
the two possible values for in Eq. ([5|). 

For Z) > 4 detectors, it is not straightforward to determine the equation of 
the admissible surface. However, we notice that any triplet taken from the set 
{Ti(0, 0)}i=2,...,_D must verify the corresponding ellipsoid equation. The admissible 
surface is thus the intersection of C^^~^ ellipsoidic cylinders. Those equations are 
redundant. As a sky location is associated to one travel time triplet only, all travel 
times can be determined from three of them. We need one equation to get one "master" 
travel time triplet and 15 — 1 — 3 equations to determine the remaining travel times, 
thus a total oi D — 2> equations. 



Sparse sky grid for the coherent detection of gravitational wave bursts 



6 



3. Sampling the travel time space and mapping to sky grid 

The coherent analysis follows directly from the computation of the "global" or 
"network" likelihood ratio testing the presence of a signal in the output from all 
detectors jointly. The network likelihood ratio involves linear combination of the data 
streams such that the response of the detectors to an incoming gravitational wave adds 
coherently. This imphes the adjustment of the data streams in time and phase [3 HI [7]. 
Time delays are applied to "synchronize" the various responses, i.e. compensate the 
travel time of the wave between detectors. Clearly, these time delays are intimately 
related to the travel time coordinates introduced in the previous section. 

Following the principle of the generalized likelihood ratio test, in presence 
of unknown parameters, the network likelihood ratio (more precisely its — log) is 
maximized over those parameters. Here, the parameters are the coordinates of the 
source location (f> and 9 and the parameters connected to the physical characteristics 
of the source (e.g., orientation of the orbital plane for an inspiralling binary, masses 
of the binary stars) and to the waveform morphology (e.g., central frequency for a 
sine-Gaussian type bursts). 

We are interested here in the maximization over cj) and 9. Mathematically, this 
is a non-linear optimization problem (due to the adjustment of the time delays in 
particular). The standard approach for its resolution is to find the maximum of the 
network likelihood ratio over a sky grid. Instead of the usual scheme based on the 
uniform discretization of the spherical coordinates, we propose here to discretize (the 
admissible surface of) the travel time coordinates introduced in the previous section. 
The resulting sky grid is built using (at least part of) the geometry of the detector 
network (i.e., the relative position of the detectors) and is thus more adapted to the 
problem at hand. The procedure goes in two steps. We first present the adequate 
sampling of the travel times. Then we map the sampled travel times to actual sky 
locations. 

3.1. Sampling the travel time admissible surface 

The data streams at the detectors are sampled at a given sampling frequency fs . We 
propose here to sample the travel time coordinates with the same pace. This presents 
the great advantage that the time delays applied in the coherent analysis are an 
integer number of samples. No interpolation procedure is thus needed to compute the 
coherently combined data streams. 

In the following, we will label with a superscript s the sampled variables. Let tf 
denote the closest time sample from time ti at detector i. We have tf = tifs+ with 
€i e (—0.5, 0.5) the truncation error. 

The travel time being the time difference — ti — ti, its sampled version 
T? = tf — if is affected by two truncation errors since r/ — Tifs + ei — ex. Now 
considering the entire set of coordinates defined in Sec. 12.11 we have 

T'^Tfs+e-ei\ (10) 

where 1 is a vector with all the entries equal to 1. 

For all T on the admissible travel time surface, we want to find the integer vectors 
provided that ei and the components of e are in (—0.5, 0.5). Let us first assume that 
ei = 0. We accept if the admissible surface intersects the cube of edge 1 centered 
on T*. Now, consider that ei ^ 0. We accept t'^ if there exists ei G (—0.5,0.5) such 
that the admissible surface translated by eil intersects the cube of edge 1 centered 
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Figure 2. Admissible surfaces of travel times for a three (left, ellipse) 
and four (right, ellipsoid) detector network and their sampling (/s = 1024 Hz). 
The unit is milliseconds. The positions of the detectors corresponds to the 
actual positions of Virgo (Pisa, Italy), LIGO H (Hanford, WA, US) and LIGO 
L (Livingstone, LA, US) with GEO (Hannover, Germany) added for the four 
detector case. The sample points obtained with the method described in Sec. 13.11 
are indicated with blue dots. Part of the points appear to fall slightly off the 
admissible surface due to round-off (see discussion in Sec. I3.2|l . 



on . This method ensures that the boundaries of the admissible surface get well 
sampled. This precaution is essential for four and more detector case since all points 
of the admissible surface belong to the boundary. 

We demonstrate this method for D — 3,4 detector network as illustrated in 
Figure [21 For a three detector network, the admissible surface is an ellipse. Most 
of the admissible samples lie inside the ellipse for which the proposed method is very 
efficient. For D=4, the admissible surface is an ellipsoid which makes the admissibility 
criterion difficult to verify. This is mainly because all of the admissible points lie on 
the boundary and the criterion amounts to determining the intersection of an ellipsoid 
and a 3-cube. However we simplify the selection criterion by looking for cubes that 
have at least one vertex inside the ellipsoid and at least one outsidqj. The result of 
this simplified method is presented in Figured 

For larger detector network, this basic method has to be abandoned since the 
admissible surface is not known analytically. One possibility is to use Monte Carlo 
(MC) trials: generate a random set of uniformly distributed sources, compute the 
wave arrival times and collect the list of distinct discretized travel times. The MC 
method has the additional advantage that the samples with a very low probability of 
appearance, i.e. samples with a small intersection between the sample cube and the 
set of translated admissible surface, are automatically dismissed. 

Table [1] gives a comparison of the different discretization methods. As expected, 
taking the physical properties of travel times (the admissible surface) into account 
decreases significantly the number of samples. Also the larger the detector network, 
the more samples are discarded by the MC method. This means that the number of 
very unlikely samples increases with the number of detectors. Note finally that the 
proportion of samples being dropped by the MC method increases with the sampling 
frequency. 

I It is possible to refine this criterion by considering more points on the surface of the cube instead 
of the vertices only 
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512 


1024 


2048 


4096 


A: 3-det 




841 


3,249 


12,543 


49,275 


Na 


285 


1,002 


3,726 


14,341 


Nmc 


285 


1,001 


3,723 


14,336 


B: 4-det 




4,205 


29,241 


188,145 


1.43 X 10'^ 


Na 


869 


4,101 


22,611 


143,155 


Nmc 


866 


3,983 


16,650 


66,454 


C: 5-det 




966,735 


1.35 X lO"" 


1.96 X 10» 


3.07 X 10" 


Nmc 


5,725 


22,985 


90,740 


351,171 



Table 1. Comparison of various for the discretization schemes used for 
the travel time space. — Nc, Na and Nmc arc the total number of grid points 
of the travel time grids obtained (i) if their physical properties are ignored (i.e., 
uniform sampling of between [— T,;, T^]), (ii) if we sample the admissible surface 
with the method described in Sec. 13. II (iii) from the Monte-Carlo search described 
in Sec. 13.11 resp. The network A is composed of Virgo, LIGO H and LIGO L. 
The network B is composed of these 3 detectors plus GEO (Hannover, Germany). 
The network is composed of these 4 detectors plus a detector located in Japan 
(assumed to be at the location of the decommissioned detector TAMA.) 

3.2. From travel time coordinates to spherical angles 

If physically admissible, travel time coordinates can be mapped into spherical 
coordinates using Eqs. (|4]) and ([5]). Some of the travel time samples fall close but 
outside the admissible surface because of the sampling round-off error (see Fig. [2]). 
While this seems marginal in the 3-detector case (because it happens to the small 
number of points at the boundary of the admissible ellipse), it is the case for almost 
all samples in the 4-detector case. A procedure is needed to map points that are 
slightly off the admissible surface to a sky position. We investigate this question here. 

Assuming a given , we search for the closest r belonging to the travel time 
admissible surface S: 

T* = argmin ||tV. - rH^, = Vsir'^fs), (11) 

where ||a;||x is the norm associated with the inner product induced by the positive 
definite matrix X. Vsi') is thus the projector onto S defined by this inner product. 

We may opt for a standard least-square minimization by setting X to identity. In a 
more accurate version, we may also model the timing errors in Eq. (jlOp as Gaussian 
random variables with zero mean and variance af and use X — diag[(T|, . . . , a^] + 
afll'^. In this last case, r* in Eq. (fTTjl is the maximum likelihood estimator of r 
(and associate sky location) . The use of this estimator does not restrict to the present 
sky sampling problem, but it is also useful in the standard triangulation scheme |llj 
used to obtain the source position from a set of measured arrival times. Note that the 
model (i.e., the values of ai) may be tuned to integrate various timing errors (such as 
calibration uncertainties). 

Let us consider networks with D = 3 or 4 detectors. In the case where X = I 
and including the quadratic expression of S in Eq. ([6]), we are led to the following 
constrained minimization problem 

minr^T - 2(tVs)^t subject to t'^At = 1 (12) 

T 

where A = Ad/Bjj. The case with an arbitrary X can be re-casted into this one by 
setting r' = X^/^r and A' = {X'^^/^)'^ AX.^/^ . This problem can be solved with the 
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method of Lagrange multipliers. It follows that the solution is r (I + A^) ^t'^ fs 
where the Lagrange multiplier A satisfies the polynomial equation: 

{T^f.fil + XA)-^AiI + XA)-'T\f, = 1. (13) 

This equation can be solved by standard root finding algorithms (in the 3-detector 
case, Eq. ([T3| is a fourth-order polynomial that can be solved analytically by radicals 
with Ferrari's method). 

This methodology has to be adapted for networks with D > 5 detectors. As 
discussed in Sec. 12. li the admissible surface is described by a set of (at most) 
D — 3 quadratic equations (each corresponding to the ellipsoid equation for a travel 
time subset). Consequently, the number of constraints imposed in the minimization 
problem of Eq. (fT^ is increased to D — 3. This problem cannot be solved using 
the Lagrangian multipliers as before. It can however be viewed as a quadratic 
programming problem with non-convex quadratic constraints for which numerical 
methods have been developed [121 [131 E] ■ 

We have applied these mapping algorithms to the travel time grids presented in 
Fig. [2] with results in Fig. [3] and to a five-detector network with result in Fig. [H 
Note that the spherical maps in these figures as well as all the other presented in this 
article are obtained using a Hammer- Aitoff projection and using as reference frame 
the geographic coordinate system. 

As expected the sampling is not uniform as it depends on the network geometry 
which is irregular. In particular, the grid gets loose in the vicinity of detector plane 
(<^ = in the reference frame in Fig.[T]). The 3-detector sky grid[^ appears to be much 
more structured than the 4-detector one. An explanation of these structures is given 
in the caption of Fig. [3l 

4. Smallest grid with complete sky coverage 

In the previous sections, we have shown how to produce a sky grid from the 
discretization of the travel time/time delay space. As this grid is built from timing 
information only, it does not use all the available information (i.e., detector orientation 
and aperture) . The consequence is that this grid is "oversampled" . Its large 
size generally prevents coherent searches to be performed in real time with current 
computers. We want to extract from the sky grid the smallest subset ensuring the 
entire sky coverage. This means that we want that the loss due to the mismatch 
between any source direction and the closest grid point be not too large. In this 
section, we describe the selection procedure of the sky grid which tightly covers the 
sky sphere with the minimal number of grid points. Our procedure results in selecting 
the grid points in areas of the sky where the position estimation is accurate, and 
withdrawing samples in areas where the accuracy is poor. 

4.1. Angular resolution from the network statistic expansion 

Coherent analysis follows directly from the computation of the network likelihood ratio 
testing the presence of a signal in the output from all detectors jointly. The likelihood 

§ For the 3-detector network, the grid points that are strictly inside the admissible ellipse have been 
mapped to the two opposite sky directions owing to the sign ambiguity in Eq. (|5]l, whereas the points 
located on and outside the boundary of the ellipse (they correspond to ^ = i.e., no sign ambiguity) 
are mapped to only one sky direction. 
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Figure 3. Sky grids associated to the sampling of the admissible travel 
time surface presented in Fig. [2] obtained by the robust mapping proeedures 
described in Sec. 13.21 The grid obtained for the 3-detector network (left) shows 
a much more structured distribution than the one obtained for the 4-detector 
network (right). This structure can be briefly explained as follows: first, the 
"sinusoidal" line which splits the figure in two halves corresponds to the detector 
plane ( </< = in the network frame, see Fig.[T]l. The two crosses X (in red) indicate 
the points {(p,0) = (0,0) and (0, tt). The grid results from the intertwinement of 
several collections of concentric circles: the first set of circles is contours of the 
polar angle (or equivalently T2, indicated in red), the second refers to contours 
of T3 (semicircles shown with black +) and the third set is formed by contours of 
T2 = Ts + C*. The two sky points along which the third set is centered correspond 
to the minimum and maximum values of C*. Geometrically, for a special case 
of VHL where T2 ~ T3 , these 2 sky points correspond to the sky directions along 
the axis (in the LHV plane) perpendicular to the angle bisector D2 — D\ — Dz 
(or L-V-H). It mimics the geometry of the sphere with these two points acting as 
poles and the contour T2 = T3 acts like an equator. In any three detector case, 
this equator is defined by = and 9 = tan~^(T2 — T3 cosq!23)/(T3 sin«23). 






Longitude 



Figure 4. Sky grid associated to the sampling of the admissible surface 

for the five-detector network C (see Table [Til and assuming fs = 1024 Hz. Clearly, 
the sky sampling is much tighter than the one obtained for other networks. This 
is due to the far location (Japan) of the additional detector. The information 
given by the time delays between this detector and the others is very accurate 
and thus valuable. 
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ratio is a function L of the data and of the unknown parameters that characterise the 
expected signal. According to the principle of the generalized likelihood ratio test , 
the network likelihood ratio (its — log, precisely) is maximized over the parameters. 
We assume here that this maximization has been done for all parameters but the 
source coordinates 4> and 9, which leads to the likelihood sky map L{(t),9). When a 
signal is present, it is expected that this function peaks where the source lies. In this 
section, we investigate the width of the peak. 

Let us consider that there is a GW source in the direction (00, ^o)- We obtain 
the peak width by the Taylor expansion of L{(j), 0) about ((/)o, 6*0): 

L(</), 9) « i(0o, ^0) + ^ ^ a2^i|^„,eo AaA/? 

a, 13 

where a, (3 — cj) or 9 in the summation and we have A(/i = cj) — (/jq and A9 — 9 ^ 9q. 

In the case where the incident wav <§ hit) = A{t) expi(p(t) is a quasi-periodic 
burst of amplitude A{t) and phase (p{t), the likelihood ratio is expressed as [5] 

D 

m9)^ J2 PnmCnC:,, (14) 
n,m— 1 

where the correlation measurement is C„ = / x„{t)h*{t — Tn)dt, denoting Xn{t) the 
strain data at detector n and h{t) the expected waveform defined above normalized 
to unit £2 norm. The operator P = {P)mn projects onto the vector space generated 
by the antenna pattern vectors F = (F)„ and F*. 

The second derivatives of the likelihood ratio can be calculatecQ following [5] 
(Sec. V A). The Hessian matrix reads 

D 
n.m — l 

where the weights are Pmn = PrnnFnF*^. The metric components can be calculated 
explicitly 

ffS" = j j A\t)A^{s)^^{^{t)Tr,-^{s)Tra)^p{^{s)Tn-mrv^)dtds{l&) 

by invoking slow variations of the amplitude as to compared to the phase's A/ A <^ ip 
and denoting A{t) the signal amplitude after whitening (i.e., division of the spectrum 
by \/ S{f) where S{f) designates the PSD of the instrumental noise which we assume 
identical at all detectors). 

If the burst is monochromatic ip{t) — 271 f^t + ipQ with a constant amplitude 
A{t) = 1 (which the case we consider in the simulations presented in the next sections), 
this simplifies to 

ff™ = ^0 [dc. ■ dp] (2^/0 (t„ - t,„ ) ) , (17) 

II We assume here that the wave is circularly polarized. However, it is straightforward to 
accommodate cases where polarization is elliptical e.g., inspiralling binary with arbitrary inclination 
w.r.t. the line of sight. 

^ For simplicity, we give the expression of Cn assuming white Gaussian noise. It is straightforward 
to extend this expression to the colored noise case by computing the correlation in the frequency 
domain and dividing by the noise power spectral density function. 

+ As in [5], we assume that the variation of the likelihood ratio is dominated by the effect of the time 
delays. We thus neglect the contribution from the phase shift due to the varying antenna pattern. 
* We used the convention ti = 0. 



Sparse sky grid for the coherent detection of gravitational wave bursts 



12 



where oc l/S{fo). It is interesting to note that the diagonal terms g^'^ are zero so 
that the Hessian matrix in Eq. vanishes when the projection operator is diagonal. 
This is true for instance when the number of detectors is D = 2 as P = I in that case 
(in words, this is equivalent to saying that it is not possible to locate the GW source 
from two detectors only). 

In the case of co-aligned (but not colocated) detectors (with identical antenna 
patterns F), the above expression are further simplified. We realize that this 
configuration is unrealistic as the detectors are built on a spherical Earth. However, 
it provides us a case where calculations and simulations are simpler. In this case, 

Pmn = \F\yD 

dl^L\Mo = -^0^(2^/0)' - ^Mirn - Tm) (18) 

nrn 

With the above Taylor expansion, we can estimate the SNR loss due to a mismatch 
between the pointing direction and the exact source direction. Let ^ denote the 
tolerable fraction of amplitude SNR loss. The grid cell is thus defined by 

J2 9l^iUo.«o AaA/3 < 2L{^o, 0o) (l - m') (19) 

a, (3 

withL{^o,Oo)^A^\F\\ 

This equation represents an ellipse centered at (0o, ^o)- The solid angle subtended 
by the ellipse expresses the angular resolutior^ achievable at this particular sky 
position. In Fig. [5l we show the (logarithm of the) angular resolution computed 
for a pure tone. It is useless to sample the celestial sphere with a bin size smaller than 
this ellipse if the goal is to get a fractional SNR loss of at most fj,. Note that while this 
requirement (with typically /i ~ 95%) is adequate for detecting a potential source, it 
is probably not sufficient for estimating its sky position. 

4-2. Smallest grid extraction as a covering problem 

We want to cover the sky sphere with the smallest number of adjacent sky resolution 
ellipses. Let He = {sk = {4>k,dk)} be the sky grid obtained from Sec. 13.21 and let £k 
be the set of samples of ile located inside the sky resolution ellipse defined in Eq. 
associated to Sk- 

We define that two grid points Sj and Sk are adjacents if Sj £ Sk and Sk G £j- 
Note that the symmetry of this criterion is important as Sj G £k does not imply that 
Sk G £j since the sky resolution ellipses change from points to points. 

We search for the smallest subset of samples ilg such that all samples in fie 
are adjacent to at least one sample in fi^. This problem is similar to the template 
placement problem encountered for the search of coalescing binaries [16l [IT] but there 
are two major differences. First, an initial sky-grid f2e is imposed here. Second, the 
topology of the parameter space to cover is different: the periodicity of the sphere 
makes the covering more difficult. In fact, there is no known optimal solution to this 
problem in the general case. Only an approximated solution can be constructed. 

(j Note that the present use of the term "angular resolution" is somewhat non-standard. Here, we 
define the width of the peak relatively to its maximum. It is therefore independent of the amplitude 
of the gravitational wave. This quantity is equal to what is referred to as "angular resolution" up to 
a scaling factor depending on the signal-to-noise ratio and fi. 
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Figure 5. Estimate of the angular resolution given by the (logjg) solid 
angle (in square degrees) subtended by the ellipse defined in Eq. (|19|l 

for the 3 detector network A in Table[T]assuming coaligned detectors (left) and real 
orientation of the detectors (right). The signal is a pure tone of frequency /o = 150 
Hz. In the first case, the angular resolution degrades when the coming wave 
has a small incidence w.r.t the detector plane (a similar eflfect appears with the 
classical triangulation method). The situation is clearly more complicated in the 
second case: the points where the resolution degrades form complicated patterns 
resulting from the interplay between the relative positions and orientations of the 
detectors. Note that the median resolution is about 9 square degrees when the 
detectors are assumed to be coaligned, while it is 10 times larger, about 94 degrees 
when considering their true orientation. 

4-3. Solving the covering problem by a greedy procedure 

The problem we have at hand can be formulated as an integer optimization problem 
with linear constraints (see [Appendix A| Eq. (|A.l|l ). In this form, it is a particular 
case of the Set Covering Problem (SCP) jTBllT^ED] which is known to be NP-complete 
|21| : no polynomial-time algorithm can solve it. Many methods have been proposed 
(for a survey, see [22]) and we selected the ones of [23] and [24] that can handle 
large number of variables (typically, the size of fig is ~ lO'* to 10^) and constraints 
(equal to the square of the initial grid size). We describe the algorithm in details in 
[Appendix A[ It consists of an efficient greedy procedure based on a dual Lagrangian 
relaxation of the problem. The greedy procedure iteratively builds a good solution: 
at each iteration, the variable Xk that maximizes a specific cost function is added to 
the solution and all rows j such that Xj and Xk are adjacents are removed from the 
problem, producing a new reduced problem that can be treated in a similar way. In 
[Appendix A[ we highlight the modifications we made to adapt to the present problem. 
In particular, we use a cost function based on the information obtained from the sole 
dual Lagrangian relaxation of the linear problem. 

4.4- Application 

In this section, we illustrate the method we propose by applying it to two different 
networks. We consider fictitious networks composed of detectors with parallel 
orientations. We can check our results against estimates of what should be the grid 
size in this case. We also consider realistic networks using the actual position and 
orientation of the main detectors currently in operation. 

From Eq. (fTS)) . we get the sky resolution ellipses for every samples assuming that 
the GW has frequency /o — 150 Hz. This frequency approximately corresponds to the 
best sensitivity for LIGO and Virgo detectors. We build the corresponding adjacency 
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fs (Hz) 


1024 


2048 


4096 


A: 3-det 


exh. set (|rie ) 


1,868 


7,184 


28,153 


min. set (|ris|) 


1,774 


1,127 


1,284 


m/\n,\ (%) 


95.0 


15.7 


4.6 


B: 4-det 


exh. set (|5^e ) 


3,983 


16,650 


66,454 


min. set (|ris|) 


1,454 


1,180 


1,354 


(%) 


36.5 


7.1 


2.0 


C: 5-det 


exh. set (|r2e ) 


22,985 


90,740 


351,171 


min. set (|fis|) 


2,695 


2,979 


3,368 


(%) 


11.7 


3.3 


0.96 



Table 2. Comparison between the sizes of initial sky grid f2e and 
minimum sky grid Qs assuming a parallel orientation for the detectors. We 
get a typical size of few thousands of vertices which is consistent with the rough 
estimate drawn from the median resolution in Fig. [5] 

matrix which feeds the greedy algorithm based on the dual Lagrangian relaxation. 

Considering that there are approximately 41,253 square degrees in the whole 
sphere. From the median angular resolution obtained in Fig. [51 we expect the full sky 
coverage with a grid of 4500 vertices in the case with coaligned detectors and about 
450 vertices considering the true orientation of the detectors. 

Detectors with parallel orientations We first examine the three-detector network A 
in Table [1] (Virgo, Ligo H and Ligo L, assuming identical noise PSDs and identical 
orientations) and four detector B (adding GEO). The samphng frequency is set to 
fs = 4096 Hz. The size of initial sky grid produced by the procedure in Sec. [3] is 
I fie I = 28153. The set covering problem can be qualified as large as the constraint 
matrix in Eq. (|A.1[) is of size \fte\ x \ile\- The application of the greedy algorithm 
based on the dual lagrangian relaxation results in the sky grid presented in Fig. [6] 
which contains \^ls\ = 1284 vertices, i.e. ~ 4.6% of the initial number of samples. 
This result is consistent with that of [5]. Rescaling the results obtained for case III 
of Table (II) to the special case discussed here, we get a sky grid of size ~ 1780 
which is comparable to jfisl mentioned above. 

This figure also displays a zoom of the resulting grid along with the corresponding 
ellipses in two small areas. In the area where the ellipse shape and orientation remain 
roughly constant, we see that the coverage is regularly performed by adjacent ellipses. 
In area where the ellipses are stretched and have orientation that changes rapidly (this 
is the case close to the line cj) = in the detector plane), some degree of overlap is 
required to have the complete coverage. This case illustrates the importance of the 
symmetric adjacency proposed and used here. 

The minimal sky grid has been computed with various detector networks and 
various sampling frequencies. The results are tabulated in Table [5J As expected the 
size of the minimum grid remains more or less constant with the sampling frequency, 
as it is prescribed by the angular resolution achievable by the detector network. 

To check that good detection performance can still be achieved with the minimum 
grid, we also compared the performance in terms of false alarm and detection 
probabilities of both the initial grid and the minimum grid. False alarm and detection 
rates were evaluated by Monte Carlo method using noise only (using simulated white 
Gaussian noise) and signal+noise trials. A monochromatic signal with frequency 
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Figure 6. Minimal sky grids for the 3 detector network A (top left, see 
Table [l] for the description of the networks) and 4 detector network B 
(top right) assuming parallel orientations (sampling frequency fs = 4096 
Hz) and signal frequency /q = 150 Hz. The number of vertices in the grid is 
1284 in the first case and 1354 in the second. Bottom row, we present two blow-up 
for the 3 detector case: (left) where the sky resolution ellipses have constant size 
and remain aligned and the other; (right) where the ellipses are strechted and 
vary rapidly from samples to samples. 



/o — 150 Hz was injected in the 3-detector network A using a sampling frequency 
fs — 4096 Hz. The position of the source was uniformly drawn over the sky sphere. 
The other signal parameters such as polarization or amplitude were also randomly 
chosen. 

We used the (-log) likelihood ratio introduced in Sec. 14. li as the detection statistics 
following the implementation given in ^. This corresponds to applying a matched 
filter filtering procedure which results in a likelihood sky map. The maximum of this 
likelihood map was compared to a threshold to decide for the detection of the GW. 

These simulations were repeated for different thresholds. For threshold values, 
10^ Monte Carlo simulations were run. The result obtained is presented in Fig. [T] We 
observe that the performances obtained with the minimal sky grid are comparable to 
that of the ( "oversampled" ) initial sky grid. 

Detectors with their true orientations Now, we consider the detectors with their 
true orientation. We show in Fig. [8] the result we get for the four-detector network 
B in Table [TJ The size of the minimum grids are gathered in Table [3] for various 
configurations. 

The resulting grids contain less samples than in the previous case (parallel 
detectors) due to the additional loss of accuracy which arises when the detectors 
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0.01 

false alarm probability 



Figure 7. Detection performance in terms of false alarm and detection 
probabilities for the 3-detector network A (sampling frequency fs = 4096). 




Figure 8. Minimal sky grid (887 
B (see Table [ij Virgo, Ligo H, 
orientation (sampling frequency 
fo = 150 Hz.) 



vertices) for the 4 detector network 
Ligo L and GEO with their actual 
fs = 4096 Hz and signal frequency 



don't share the same orientation. 



5. Conclusions 

Most coherent detection scheme includes an algorithmic loop over source directions 
taken from a predefined sky grid. This part of the detection algorithm takes a 
significant amount of computing ressources. One of the available degrees of freedom to 
address this issue is the way the sky sphere is sampled. We investigated the possibility 
to build an optimally "sparse" sky grid which performs the complete sky coverage with 
the smallest number of vertices. It is clear that the optimal sky grid should adapt the 
pointing accuracy. 

Contrarily to electromagnetic antenna arrays (e.g., used for RADAR applica- 
tions), gravitational wave detector networks are not arrays of regularly spaced and 
oriented sensors. The pointing accuracy we get with such network results from the 
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fs (Hz) 


1024 


2048 


4096 


A: 3-det 


exh. set (|rie ) 


1,868 


7,184 


28,153 


min. set (|ils|) 


429 


391 


440 


(in%) 


23.0 


5.4 


1.6 


B: 4-det 


exh. set (|5^e ) 


3,983 


16,650 


66,454 


min. set (|ris|) 


937 


817 


887 


\ns\/\n,\ (%) 


23.5 


4.9 


1.3 


C: 5-det 


exh. set (|5^e ) 


22,985 


90,740 


351,171 


min. set (|fis|) 


1,587 


1,723 


1,898 


(%) 


6.9 


1.9 


0.54 



Table 3. Comparison between the sizes of initial sky grid f2e and 
minimum sky grid Sis considering the actual orientation of the detectors. We 
get a typical grid size of few hundrends of vertices which is consistent with the 
estimate drawn from the median resolution in Fig. [5] 

complicated interplay between the relative positions and orientations of the detectors. 
Because of the network heterogeneity, the pointing accuracy varies largely over the 
whole sky. In this paper, we build sky grids from an estimate of the "local" angular 
resolution. The vertex density is reduced where the resolution is coarse and vice- versa. 

The method goes into two steps. We first produce a grid from pre-selected sky 
directions where the travel time between detectors and the time reference is (close to) 
an integer number of time samples. Those points are convenient because the time- 
shifts necessary to align in time the signals (i.e., compensate the travel time between 
the detector and the time reference) received by each detector are trivially performed. 
This operation amounts to sampling the manifold described by the travel times. This 
sampling can be determined analytically in the case of three (i.e., the manifold is an 
ellipse interior) and four detectors (i.e., the manifold is the surface of an ellipsoid). 
For larger networks, we suggest to apply Monte Carlo procedures. Once done, a key 
point is to map from travel times to sky directions. To perform this mapping, we 
propose a method that is robust to the round-off error due to the truncation of the 
travel times. 

In the second step, we extract from the first grid the smallest sub-set of vertices 
that ensures complete sky coverage for a given loss in signal-to-noise ratio. We build 
this sub-set by casting the grid size minimization into a set covering problem. The 
procedure we propose to solve the corresponding linear program with linear and 
integer constraints relies on a greedy algorithm and a dual Lagrangian relaxation. 
The resulting grid show considerable reduction in size. We have checked that the 
results we get with this grid are comparable to what we get when using more resolved 
sky grids (i.e., at a given false alarm rate, the detection probability are comparable 
when performing an all-sky blind search). Our investigations have been limited to 
source detection. It is likely that the sky grid has to be refined if the goal is the 
accurate determination of the source position. 

The overall procedure to build the grid can be quite computationally heavy. 
However this procedure needs to be done only once beforehand. 
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Appendix A. Solving the set covering problem 

We propose to formulate the problem of building the minimum sample set as an integer 
optimization problem with linear constraints in the following way: 



Here the vector x of size |f2e| is the output vector that should define the minimal set: its 
entry Xk takes value 1 if the corresponding sample Sk is included in fig, and otherwise. 
Hence, the non-zero entries of x define the subset fig, and the objective function of 
this optimization problem is the number of samples in the final set (the number of 
non-zero values of vector x, to be minimized). The vector u is a column vector of size 
I fie I such that u = [1, . . . , 1]'^. Finally the matrix A of size (|rie| x l^^el) represents 
the neighbouring relationships: it is such that Qij = 1 if and sj are neighboors and 
Uij = otherwise. Each row j of matrix A represent a specific constraint on sample 
Sj in the oversampled set: it is here to insure that Sj is be adjacent to at least one 
non zero entry of x, and therefore covered by at least one resolution cell. 

As the matrix of linear constraints A may contain several tens or hundreds of 
thousands of columns and the same amount of rows (see the number of samples from 
Table [TJ the problem is very difficult and cannot be solved by optimal methods. We 
therefore propose to use here a method inspired by [53] and [23] . It follows an efficient 
greedy procedure based on a dual Lagrangian relaxation of the problem. The greedy 
procedure iteratively builds a good solution in the same way as in [24j : at each 
iteration, the variable x^. that maximizes a specific cost function is added to the 
solution and all rows j such that xj G £k are removed from the problem (because the 
corresponding constraints are fulfilled by the solution), producing a reduced problem 
that can be treated in a similar way. Here the cost function is computed thanks to 
the information obtained from the dual Lagrangian relaxation of the reduced problem. 
Contrary to the method proposed in [24] , we do not use the primal relaxation of (jA.ip , 
which reduces the computational cost of the algorithm as well as the dependency of 
the method over tunable parameters, and leads to smaller running time and slightly 
better solutions. 

The dual problem of (|A.1[) is defined as: 

max u^y 

s.t. A^y < u„, (A.2) 

yfc >0, k = i,...,\n,\. 

The Lagrangian relaxation of the dual problem (|A.2p is given here by: 



where /x is the vector of Lagrange multipliers and L(y,/x) = u^jY — /i'^(A'^y — u„) 
is the Lagrange function of the dual problem. Here the maximization of L(y,/x) over 
G {0, 1}, fc = 1, . . . , n is straightforward: indeed it is obtained by setting ?/fc = 1 
if the k^^ column of the row vector — ^^A^ is positive, and i/k = otherwise. 




(A.l) 



k=l 

S.t. Ax > u, 

Xk e {0,1}, k = l,...,\n,\. 



min maxL(y./x) 
p y 

s.t. Mi > 0, i = 1, . . . , m. 



(A.3) 



< yfe < 1, /c = 1, . . . ,n. 
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Besides the function maxy _L(y,/i) is convex in /x, and therefore it can be minimized 
efficiently by means of a subgradient algorithm. Two main problems arise: 

• the optimal values of (IA.2|) and of its Lagrangian relaxation (|A.3|) may not be 
equal (duality gap); 

• the vector y*(/x) minimizing _L(y,/x) is generally not feasible, i.e. it doesn't verify 
the linear constraints A-^y*(/x) < u„. 

However, despite these problems, and contrary to the multipliers obtained from 
the primal Lagrangian relaxation 123], the dual Lagrange multipliers /i carry very 
useful information about problem (jA.ip . In particular /x is a near-optimal and near- 
feasible solution of the linear relaxation of (|A.1[) . It can therefore be used to decide 
which variable should be added to the final solution. Here the variable chosen is 
simply the one corresponding to the largest entry of fj.. It should be also noted that 
all samples Sk such that = {s^} must be in the final solution set since these points 
have no adjacent point. 

The final greedy algorithm can be summarized as follows: 

(i) Find the oversampled set fie', 

(ii) Build the adjacency matrix A; 
(in) Set the solution set il^ = 0; 

(iv) Find all points Sk such that £k = {sfc}, add them to f2s, remove them from the 
set fie and remove the corresponding rows and columns from the matrix A; 

(v) While Qe ^ 0: 

• Solve the dual Lagrangian relaxation (jA.Sp of problem (jA.lD : 

• Find the largest Lagrange multiplier /Zfc , add the corresponding sample Sk to 
Qs, remove the corresponding column from matrix A as well as all rows j 
such that Sj and Sk are neighboors, and reduce the set fig similarly. 

It is also sometimes possible to improve the quality of the solution by introducing 
a randomized criterion for the variable selection. In [24j , it was proposed to run several 
time the algorithm. For the first run, the variable selected is, as described above, the 
one corresponding to the largest entry of (j,. However for the other runs, the variable 
selected may be the one corresponding to the second largest entry of /i with probability 
0.5. We adopted this strategy since the choice of the minimum sample set has to be 
done only once in our case and therefore the additional computational cost induced 
by running several time the same algorithm is not prohibitive. 
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